close all;clear
%给定参数
Nt = 80;
p1 = [100, 100, 10, 10, 10, 10, 100, 100];
p2 = [10, 10, 100, 100, 50, 50, 50, 50];
p3 = [50, 50, 50, 50, 100, 100, 10, 10];
h1 = 1000;
h2 = [200, 5000, 200, 5000, 200, 5000, 200, 5000];
u = 4*pi*10^(-7);
i = sqrt(-1);
%进行计算
for ix = 1:8
    for  iy = 1: Nt
        w(ix,iy) = 10 ^ ((iy+10) / 10) / h1;
        lanbta(ix,iy) =  2*pi*sqrt(2) ./ sqrt(w(ix,iy)*u./p1(ix));
        k1 = sqrt(-i * w(ix,iy) * u / p1(ix));
        k2 = sqrt(-i * w(ix,iy) * u / p2(ix));
        k3 = sqrt(-i * w(ix,iy) * u / p3(ix));
        R2=tanh((h2(ix) .* k2 + atanh(k2 ./ k3)));
        R1=tanh( h1 * k1 + atanh((k1 ./ k2) .* R2));
        P(ix,iy) = R1 .^ 2;
        A(ix,iy) = abs(P(ix,iy));
        B(ix,iy) = angle(P(ix,iy));
    end

end
%进行画图
figure(1)
loglog(lanbta(1,:)/h1,A(1,:))
title("H型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(2)
loglog(lanbta(1,:)/h1,A(2,:))
title("H型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(3)
loglog(lanbta(1,:)/h1,A(3,:))
title("K型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(4)
loglog(lanbta(1,:)/h1,A(4,:))
title("K型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(5)
loglog(lanbta(1,:)/h1,A(5,:))
title("A型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(6)
loglog(lanbta(1,:)/h1,A(6,:))
title("A型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(7)
loglog(lanbta(1,:)/h1,A(7,:))
title("Q型三层大地电磁测深响应曲线(中间薄)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

figure(8)
loglog(lanbta(1,:)/h1,A(8,:))
title("Q型三层大地电磁测深响应曲线(中间厚)")
xlabel("lanbta1/h1"); ylabel("pw/p1")

